/*  WECC Intrastate Reallocation with bootstrapped damages.

NOTES: All directories and paths should be set in the "Create_TableA5.do" file.

*/

clear all


import delimited `"${data_path}/Generation_and_damages_withCO2_2019-11-24_eGrid_SR.csv"', case(preserve) 

merge 1:1 zip using `"${data_intermediate}/WECCBootsFinal11232019_wide.dta"'
keep if _merge==3
drop _merge
rename zip zipcode
rename STATE state
*drop num_systems

*merge 1:1 zipcode using NumberSystems
*Note: We have 1478 zipcodes with solar panels but no damage calculations.
*So I assign damages to these zipcodes based on the damages from the previous zipcode (when zipcodes are listed in numerical order)
*replace total_damages = total_damages[_n-1] if _merge==2
*drop _merge

merge 1:1 zipcode using `"${data_path}/HousingUnits"'
keep if _merge==3
*need to drop zipcodes with missing housing units or zip codes that have 0 housing units
*Note: There are 121 zip codes with a positive number of installations but 0 housing units
*This can be observed by uncommenting the next line.
*count if single_units ==0 & num_systems>0 & !missing(num_systems)
drop if missing(single_units) | single_units==0
drop _merge

foreach var of varlist total_damages1b1-total_damages8b99 {
	gen total_damages_all = `var'*num_systems
	egen SumOld`var'=sum(total_damages_all)
	drop total_damages_all

}
gen total_damages_all = total_damages*num_systems
egen sumOldPtEst = sum(total_damages_all)
drop total_damages_all

egen old_p25 = rowpctile( SumOldtotal_damages1b1- SumOldtotal_damages8b99), p(2.5)
egen old_p975 = rowpctile( SumOldtotal_damages1b1- SumOldtotal_damages8b99), p(97.5)
egen old_mean = rowmean( SumOldtotal_damages1b1- SumOldtotal_damages8b99)

***PAUSE HERE to record estimates and CIs of avoided damages from currently installed capacity

*drop SumOldtotal_damages1b1- SumOldtotal_damages8b99
*drop old_p25 old_p975 old_mean sumOldPtEst

***Get new damages based on point estimate

*generate total damages (damages per 4kW installation * number of installations)
gen total_damages_all = total_damages*num_systems

sort total_damages single_units
egen rank=rank(-total_damages)
sort rank single_units
drop rank
gen rank=_n /*ranks zipcodes by total_damages and number of houses*/

sort state rank
bysort state: egen rankbystate=rank(rank)
bysort state (rankbystate): egen sumsystemsState = sum(num_systems)
egen stategroup = group(state)
gen id = _n /*Create new overall ranking for purpose of calling values later*/

* Simulation 1 (moving all panels to highest total_damages zipcode in state subject to #of housing units)
*drop remaining
gen state_new_num_systems=0
gen remaining = sumsystemsState

sum stategroup, meanonly
local numstates=r(max)
forvalues i=1/`numstates' {
	sum rankbystate if stategroup==`i', meanonly
	local maxRankInState=r(max)
	forvalues j = 1/`maxRankInState' {
		sum id if stategroup==`i' & rankbystate==`j', meanonly
		local index1= int(r(min))
		quietly replace state_new_num_systems=min(remaining[`index1'],single_units[`index1']) if id==`index1' & remaining[`index1']>0
		quietly drop remaining
		quietly bysort state (rankbystate): egen sumNewByState = sum(state_new_num_systems)
		quietly gen remaining=sumsystemsState-sumNewByState
		quietly drop sumNewByState
	}
}

gen newState_total_damages=state_new_num_systems*total_damages
* check
bysort state (rankbystate): egen sumNewByState = sum(state_new_num_systems)
egen sumNew=sum(state_new_num_systems)
egen SumNewDamagesState=sum(newState_total_damages)

***BOOTSTRAP Damages
foreach var of varlist total_damages1b1-total_damages8b99 {
	gen boot_total_damages_run = `var'*state_new_num_systems
	egen SumBoot`var' = sum(boot_total_damages_run)
	gen gain_`var' = SumBoot`var' - SumOld`var'
	drop SumBoot`var' SumOld`var'
	drop boot_total_damages_run
}
egen gain_p25 = rowpctile( gain_total_damages1b1- gain_total_damages8b99), p(2.5)
egen gain_p975 = rowpctile( gain_total_damages1b1- gain_total_damages8b99), p(97.5)
egen gain_mean = rowmean( gain_total_damages1b1- gain_total_damages8b99)
egen gain_sd = rowsd(gain_total_damages1b1- gain_total_damages8b99)
gen gain_pt_est = SumNewDamagesState - sumOldPtEst
drop gain_total_damages1b1- gain_total_damages8b99

putexcel set `"${data_final}/tableA5.xlsx"', sheet(sheet1) modify

putexcel B3 = gain_mean/1000000
putexcel C3 = gain_sd/1000000
putexcel D3 = gain_p25/1000000
putexcel E3 = gain_p975/1000000

